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Charge hydration asymmetry (CHA) - a characteristic dependence of hydration free energy on the 
sign of the solute charge - quantifies the asymmetric response of water to electric field at microscopic 
level. Accurate estimates of CHA are critical for understanding of hydration effects ubiquitous in 
chemistry and biology. However, measuring hydration energies of charged species is fraught with 
significant difficulties, which lead to unacceptably large (up to 300 %) variation in the available 
estimates of the CHA effect. We circumvent these difficulties by developing a framework which 
allows us to extract and accurately estimate the intrinsic propensity of water to exhibit CHA from 
accurate experimental hydration free energies of neutral polar molecules. Specifically, from a set of 
504 small molecules we identify two pairs that are analogous, with respect to CHA, to the K + /F“ 
pair - a classical probe for the effect. We use these “CHA-conjugate” molecule pairs to quantify 
the intrinsic charge-asymmetric response of water to the microscopic charge perturbations: the 
asymmetry of the response is strong, ~ 50% of the average hydration free energy of these molecules. 
The ability of widely used classical water models to predict hydration energies of small molecules 
correlates with their ability to predict CHA. 


I. INTRODUCTION 


The water molecule, H 2 O, is among the simplest chem¬ 
ical structures, yet the liquid phase of water is among 
the most complex liquids. Many of its unique, vital for 
Life properties, including the ability to establish com¬ 
plex hydrogen bonded structure [ 2 ], are due to the com¬ 
plexity of electrostatic interactions. The ability of water 
to hydrate ions and biomolecules is crucial to biologi¬ 
cal functions; detailed atomistic understanding of these 
functions is impossible without accurate description of 
the energetics of aqueous solvation (hydration). Given 
the very polar nature of water as a solvent, understand¬ 
ing key details of molecular hydration requires accurate 
and comprehensive experimental characterization of elec¬ 
trostatic properties of water molecule in liquid phase at 
microscopic level. This characterization is currently lack¬ 
ing: for example, even the value of the dipole moment of 
water molecule in liquid phase is still debated the ex¬ 
perimental range is from 2.7 to 3.2D, and higher 

multipole moments are not available from experiment at 
all. 

Experiments further lack consensus in quantifying 
the key characteristic of the asymmetry of water re¬ 
sponse to microscopic fields of solvated ch arges - the so 
called charge hydration asymmetry (CHA) |6H20j. Fig- CD 
One of its earliest (21.] known manifestations is the ob¬ 
served differences in hydration free energies of oppositely 
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charged ions of similar size, e.< 7 ., the K + /F _ pair, Fig. [I] 
completely unaccounted for in the linear response con¬ 
tinuum framework (22|- These CHA-related differences 
are larger than many relevant biomolecular energy scales, 
such as folding free energy of a typical protein. The CHA 
effects are strong not only for the charged, but also for net 
neutral solvated structures [3,13, HQj • Absolute single- 
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FIG. 1: Uncertainty in experimentally determined charge 
hydration asymmetry (CHA) for K + /F _ ion pair. Shown 
are the hydration free energies, AG, for K + (top mark¬ 
ers) and F _ (bottom markers) from four classical experimen¬ 
tal references [3, I3HE], from historical to modern, span¬ 
ning more than eight decades. Also shown are the abso¬ 
lute CHA, A AG = AGTK+) - AG(F“), and the dimen¬ 
sionless relative CHA |l6j], which quantifies the asymme¬ 
try of the charge hydration: rf = | |, where (AG) = 

I (AG(K+) + AG(F-)). 

ion solvation free energies are conventionally deduced 
from experimental estimates of hydration enthalpies of 
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ionic solutions via the Born-Haber cycle [2J, [23 using 
highly controversial estimation of absolute proton (H + ) 
hydration enthalpy. Uncertainties of up to 16 kcal/mol 
are reported that stem from various degrees of extra¬ 


thermodynamic assumptions [23|, [25|, [27] or cluster ion 
solvation approximation [24] made in obtaining these es¬ 
timates. Furthermore, ion hydration free energy is di¬ 
rectly affected by liquid-vapor and liquid-cavity surface 
(interface) potentials: every time a charge q crosses such 
interface, the electrostatic energy of the charge changes 
by q&nt > where <pi nt is the electrostatic potential differ¬ 
ence associated with a the interface. However, available 
experimental estimates of this elusive quantity can differ 
by a factor of five [28], [29| ; there are fundamental difficul¬ 
ties associated with measuring and even interpreting [30j] 
the water surface potential [3l], |33 ■ 


The net result of these, and possibly other problems, 
is an almost 300% variation between the four available 
comprehensive experimental data sets of ion hydration 
energies in the strength of the CHA effect for K + /F“ 
pair, Fig. [Q At the same time, even uncertainties of 
the order of just several fcsT can have profound effect 
on the thermodynamics or kinetics of different phenom¬ 
ena where ion dehydration is involved, e.g., membrane 
ion channel transport, ion binding to receptor proteins, 
etc. An accurate quantification of the inherent CHA of 
water is therefore of paramount importance for accurate 
quantitative estimates in biology, chemistry, and physics. 


In contrast to ions, direct measurements of experimen¬ 
tal hydration energies [33l, ‘34] of neutral solutes are rel¬ 
atively straightforward and very accurate, with errors 
within a fraction of kcal/mol. The two primary sources 
of error in the case of ionic solutes are completely elimi¬ 
nated in the case of net neutral solutes, qt = 0. First, 
the contribution of the interface potential to the energy of 
the solute transfer is zero, (j>i n t ]Tb qi = 0. Second, exper¬ 
imental hydration energy estimates of neutral molecules 
are typically obtained using calorimetric measurements 
of transfer free energy, which are free from empirical ref¬ 
erence energies such as the proton hydration free energy 
needed in the case of ions. j35[ Connection to computa¬ 
tional models is also easy to make since the relationship 
between measured and computed hydration energies is 
straightforward for neutral solutes HI, l36l - l39j . in con¬ 
trast to ions, see e.g. Ref. [38] for a detailed overview 
of the associated issues. We therefore propose to quan¬ 
tify the CHA effect by utilizing experimental hydration 
energies of small neutral molecules, rather than ions. 
Specifically, we identify “CHA-conjugate” pairs of anion¬ 
like and cation-like neutral polar molecules that behave 
just like the K + /F _ pair with respect to charge hydra¬ 
tion asymmetry. To characterize the asymmetry of wa¬ 
ter response to microscopic charge independently of the 
strength of the charge probe (absolute hydration energies 
of small neutral molecules are much smaller than those 
of ions) in what follows we will quantify the strength of 


the CHA effect by a dimensionless quantity rf = 


AAG 

jAGj 


where (A G) is the average hydration energy of the pair 
of solutes, Fig. [Q 


II. RESULTS AND DISCUSSION 

Without loss of generality, the hydration free en¬ 
ergy, AG, of a molecule can be decomposed into 
symmetric, AG sym , and asymmetric, A G asyrn : parts 
with respect to inversion of solute charges: AG = 
A G s v m -f a Q as ym. suc j 1 that upon sign inversion of 
the partial atomic charges AG sym remains invariant, 
whereas A G asyrn —t —A G asym . Physically, the pro¬ 
posed symmetric-asymmetric decomposition of the solva¬ 
tion free energy reflects the inherently charge-asymmetric 
microscopic response of water to perturbing electric field 
due to a solvated charge. Here, AG asym part accu¬ 
mulates all the deviations from the symmetric response 
AG syrn that is invariant upon inversion of the perturbing 
electric field; the response due to dipole polarization of 
water contributes primarily to the polar part of AG sym . 
The non-polar component of AG for molecules is also 
mostly symmetric [14| . and is similar in magnitude for 
molecules of similar size. The magnitude of AG asym rel¬ 
ative to A G sym characterizes charge asymmetry of the 
solvent response: e.g. in a hypothetical water consisting 
of molecules with a perfectly tetrahedral charge distribu¬ 
tion this asymmetric response is zero[l6|. 

For an arbitrary pair (A,B) of molecules, the exper¬ 
imental difference in their hydration energies AAG = 
A G(A) — AG(B) would contain a mix of CHA-related 
and unrelated components. Here we propose to cleanly 
isolate the asymmetric part by identifying special pairs 
of “CHA-conjugate” neutral molecules - molecules with 
same AG sym , but equal and opposite AG asym ; for these 
special pairs AG sym would cancel out in the differ¬ 
ence, while |AG asym | would combine, so that AAG = 
A G asym (A) - AG asym (B) would contain only the CHA 
effect we seek to quantify. 

Since no known experimental method can perform the 
decomposition of the total solvation energy into charge- 
symmetric and asymmetric parts, we must resort to the¬ 
oretical water models to identify such special pairs of 
molecules. However, once such a pair is identified, the 
values of A G sym and A G asym for each molecule in the 
A/B pair can be extracted from the experimental hy¬ 
dration energies; A G sym = ±(A G{A) + AG{B)) and 
\AG as ym\ = i|AG(A) - AG(B)\. The implications of 
the decomposition for development of simplified water 
models will be discussed below. 

Within the widely used family of TIPnP water mod¬ 
els 53, the propensity to cause CHA increases [H, 
E3 progressively from TIP5P to TIP3P[4l| to TIP4P- 
Ew[42|- In contrast, one of the earliest water models, 
BNS[43, does not exhibit any CHA: due to its per¬ 
fectly tetrahedral geometry, charge inversion within this 
molecule can be emulated by a set of rotations around its 
oxygen center - a sufficient condition for the absence of 
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Perfect TIP5P TIP3P TIP4P-Ew 

Tetrahedron TIP5P-E SPC/SPC-E 



FIG. 2: A family of rigid n- point water models arranged in ascending order with respect to deviation from the perfect tetrahedral 
symmetry (BNS model, zero CHA). The ability of the water models to cause CHA is proportional to the magnitude of parameter 
d that progressively breaks the specific charge inversion symmetry in this family of models. 


CHA [ 16]. This high degree of charge inversion symme¬ 
try is progressively broken in other water models shown 
in Fig. [2] by an order parameter S that has an intuitive 
geometric interpretation. A relation between S (geome¬ 
try) and CHA (energy) is established in SI, where it is 
shown that S ~ fl 2 - the cubic octupole moment of water 
molecule earlier recognized [16j, |44| - |46|| as key for breaking 
the charge hydration symmetry. This multipole moment 
is defined as SI 2 = (1/2 )(fl yyz — ^xxz) where are the 
components of the traceless octupole tensor of a water 
molecule model in the Cartesian frame with origin at the 
water oxygen center and OH bonds in the yz plane with 
z-axis bisecting the H-O-H bond angle. 

The asymmetric part of the solvation free energy, 

A G asy m can be treated as a first order perturbation (with 
respect to the symmetry breaking variable S) to the com¬ 
pletely charge-symmetric part, AG sym = A G(6 = 0); 
i.e. A G{8) ~ AG (6 = 0) + For small 5, 

Ai Qasym ^ ttAAfi oc S. To illustrate the idea within a con¬ 
crete theoretical framework we follow an earlier work (l6j 
where CHA is introduced explicitly into the Born equa¬ 
tion [22I ] : 

aG " - 0 - 7) 0 - ■ 

( 1 ) 

Here, e is the dielectric constant of water, R w is the radius 
of water molecule, q and R are the ion charge and ionic 
radius, respectively, and R s is a constant shift to the 
dielectric boundary M- Neglecting a minor difference in 
the ionic radii of K + and F _ ions, [T] results in the same 
A G sym = AG(5 = 0) for both ions and a symmetric gap 
in hydration free energies 

Aao < i ) = 2 ' ao “"< s >' = (‘ - 7 ) TrTr7) (JTTbj ’ 

( 2 ) 

which increases monotonically with 6, Fig. [3] In other 
words, for a cation/anion pair of the same size, [l] 
describes a “CHA-conjugate” ion pair with identical 
A G sym and equal but opposite AG asym . 

We propose to use these two distinct characteristics 
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FIG. 3: K + (blue dashed line) and F _ (red dashed line) ion 
pair defines the characteristic behavior of CHA-conjugate pair 
of solutes - a symmetric, monotonically increasing difference 
in the electrostatic hydration free energies with respect to 
increasing degree of the perturbation S that breaks the charge- 
inversion symmetry of the water models. Hydration energies 
are computed via “CHA-aware” Born equationjlf|, [I] with 
Rs = 0.52 A, R w = 1.4 A, and ionic radii from Ref. H^I 


of the AG(5) for a pair of a “CHA-conjugate” solutes 
as a signature feature to identify CHA-conjugate pairs 
among neutral molecules. That is, we seek pairs of small 
polar neutral molecules for which AG(<5) behaves similar 
to that of the K + /F~ pair in Fig. [3j despite differences 
in scales of hydration free energies - hydration free en¬ 
ergies of small neutral molecules are an order of magni¬ 
tude less than that of the ionic pair. In hydrated polar 
molecules, for oppositely partially charged atoms one can 
expect local cation-like or anion-like response of water. 
The contributions to the overall CHA from the hydration 
of these atoms, especially atoms of the polar groups of a 
molecule, may be quite dissimilar. In general, these con¬ 
tributions will not cancel each other completely, resulting 
in a net CHA effect. In addition to that, in some pairs 
of molecules, these contributions may combine in a way 
to produce the AG(5) pattern we seek, Fig. [3] Note that 
the very existence of such pairs is not guaranteed a priori : 
hydration energy of a molecule, as well as its symmetric 
and asymmetric parts, are influenced by a complex inter¬ 
play of screened interactions between its partial charges. 
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To illustrate the above points mathematically we use the 
multi-atomic analog of [1] derived directly from the so- 
called CHA-GB model [181]. [5liri Computational Methods. 
Within this formalism, 

A G(6) ~ -\ (l - i) Y, jSi [1 - ^ > ( 3 ) 

i,j ' l J 

( _ ^ \" 1/2 

where jjp w = I jA + RiRje 4Kii b 1 is the 

purely charge-symmetric Green function of the 
generalized Born formula E! [4f| , and = 

r r ii / r \ r jl \ ~ 

R i R^~ T ^ K i ( sgn{-£, kqk e^-Ktk] sgn^qie ^ ] 

f&B* I R,-Rs + R lL , + Rj-R a +R w 

Note that [3] has the same structure as [T| with respect 
to variation in 5, which further rationalizes the use 
of our ion-inspired criteria, Fig. El for identifying 
CBA-conjugate pairs among molecules. Bowever, in 
contrast to the case of a pair of oppositely charged 
monovalent ions of the same size, for most pairs of 
neutral molecules their |AG asym | are not equal to each 
other even if their AG sym are the same. This is due 
to a complex interplay between the self and the 

cross (charge-charge interactions, T l3 ) contributions to 
CBA within each molecules. In fact, as we show later, 
we found only two pairs of CBA-conjugate molecules 
from an available comprehensive diverse set of 504 
small, neutral molecules jipif with known experimental 
hydration free energies (see Computational Methods). 
Note that these experiments were conducted under 
conditions to ensure jH] that these molecular solutes 
were indeed neutral, and therefore avoid the issues 
related to uncertainties in hydration energies of charged 
solutes discussed earlier. The fast analytical CBA-GB 
model E! was used in the initial screening for the 
signature AAG(S) gap, and the best candidates were 
then confirmed by careful free energy perturbation 
(FEP) calculations |34| using the three TIPnP explicit 
water models in Fig. 0 see Computational Methods. 

In the search for CBA-conjugate pairs of small 
molecules, we utilize rigid, non-polarizable water models 
of the same geometry (TIPnP family) that have roughly 
the same dipole moment (less than 0.06 D difference be¬ 
tween the models) and exhibit a well studied behavior 
with respect to CBA EM!. It was shown that the dif¬ 
ferent CBA effects these models can produce are propor¬ 
tional to the values of their cubic octupole moments B 2 
and thus to the different values of the symmetry breaking 
parameter <5 inherent to these models. It is also critical 
that the monotonic increase of the AAG(<5) gap with <5 
is symmetrical, that is AG sym is practically the same 
within this family of water models. These two proper¬ 
ties of TIPnP family are essential for the identification 
of CBA-conjugate pairs of solutes. 

Our search has yielded two CBA-conjugate pairs of 
neutral, but highly polar molecules, Fig. [I] which by 


our quantitative criteria behave just like the K + /F - ion 
pair with respect to CBA, compare Fig. [I] and Fig. [3] 
The qualitative rationalization for the unique behav¬ 
ior of these two pairs with respect to CBA is illus¬ 
trated for one of them in Fig. [5] The anion-like (cation¬ 
like) molecule has one relatively highly charged, “CBA- 
dominant” atom, e.g. the amide oxygen of N-methyl 
acetamide (carboxylic hydrogen of pentanoic acid) that 
behaves like F~ (or K + ) with respect to the structuring 
of water around it. The effect of the opposite charges 
within each molecule is much more diffuse, and can not 
cancel the dominant asymmetric contribution from this 
one atom. 
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FIG. 4: Hydration free energies of two pairs of molecules 
analogous to K + /F _ pair with respect to CHA; compare with 
Fig-El The energies are computed using the CHA-GB model 
(dashed lines) [18(], and free energy perturbation (FEP) cal¬ 
culations (solid markers), see Computational Methods. Sta¬ 
tistical error bars (not shown) for the FEP calculations are 
smaller than the symbol size 


We quantify the charge hydration asymmetry of 
these two pairs by their experimental relative CHA, 
The experimental 77 * for the two CHA- 


V = 


AA G 


(AG) 


conjugate pairs of neutral small molecules (3-methyl bu¬ 
tanoic acid, N-methylacetamide) and (pentanoic acid, N- 
methylacetamide), respectively, are presented in Tabled] 


TABLE I: Relative CHA, g *, of two CHA-conjugate pairs of 
neutral molecules from the experimental hydration free ener¬ 
gies 



The physical meaning of the large, almost 0.5, ex¬ 
perimental value of 77 * is that the intrinsic charge- 
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FIG. 5: Similarity between water-structuring effects of ions 
(left) and one pair of CHA-conjugate neutral polar molecules 
(right). Shown are water-oxygen (magenta) and water- 
hydrogen (cyan) occupancy maps computed using TIP3P wa¬ 
ter model. High probability isomaps are shown using solid 
surfaces and relatively lower probability isomaps are shown as 
translucent surfaces. Atom that contributes predominantly to 
the CHA effect in each molecule is pointed out by green circle. 
The other cation-like molecule, 3-methyl butanoic acid, is an 
isomer of pentanoic acid and has similar occupancy isomap, 
see SI. Using TIP4P-Ew instead of TIP3P results in similar 
occupancy maps (not shown). 


asymmetric response of water to microscopic charges 
(atomic size CHA perturbations) is strong, comparable 
to the charge hydration energy itself. For example, ex¬ 
perimental hydration energies of pentanoic acid and N- 
methylacetamide are -6.16 and - 10.00 kcal/mol respec¬ 
tively, [33| i.e. the |(AG)| = 8.08 kcal/mol while the 
AA G= 3.84 kcal/mol for the pair. In fact, the asym¬ 
metry of the individual perturbation caused by a single 
“CHA-dominant” atom is even stronger, as the net CHA 
effect in a pair of molecules measured by rj* can be at¬ 
tenuated somewhat by the opposing contributions from 
other charges, see Fig. [5] and SI. We again stress that 
the rj* from Table Q] are free of the issues related to the 
fundamental or technical difficulties [38| associated with 
the ion hydration asymmetry: 77 * inferred from CHA- 
conjugate pairs of neutral polar molecules directly and 
accurately characterize the intrinsic asymmetry of the 
water response to microscopic electric field. As expected, 
the uncertainty in the experimental rf is small, about 2 
%, as seen from the difference between the two CHA- 
conjugate pairs. This high accuracy is in contrast to the 
almost 300 % difference between available experimental 
hydration energies for K + /F _ pair, Fig. [T] 

Looking through the prism of CHA at the develop¬ 
ment of classical fixed-charge water models during the 
past several decades, we note that, geometrically, the 
largest variation between the popular water models ap¬ 
pears along the single coordinate: one that breaks the 


perfect charge-inversion symmetry of a perfect tetrahe¬ 
dral charge arrangement assumed by early BNS model, 
Fig. [2] Judging by the ability of each water model in 
Table U to predict the experimental 77 *, we conclude 
that to the extent that the charge distribution of a real 
water molecule in liquid water can be approximated by 
three point charges, the distribution mi ght be close to 
the recently proposed water model OPC [5(J,|5Jj, Fig. [ 6 ] 
This new water model is, perhaps, the very first at- 


TABLE II: Relative CHA, 77 *, of two CHA-conjugate pairs of 
neutral molecules from experiment and explicit water simula¬ 
tions. Water models used for comparison are TIP5P, TIP3P, 
TIP4P-Ew and a recent 4-point model OPC (5Ci|, [5l| 


Molecule Pair 

Exp 


Water models 




TIP5P 

TIP3P 

TIP4P-Ew 

OPC 

3-met^it^rioic ^kl W 

*' ^ N methyl acetamiHe 

0.49 

0.10 

0.43 

0.56 

0.45 

pentanoic acid V 

N methyl acetamioe 

0.48 

0.11 

0.42 

0.53 

0.45 


tempt to explicitly consider the CHA effect in water 
model optimization [5l|; the model is curiously success¬ 
ful in accurately reproducing most properties of liquid 
water [5l| . Importantly, the accuracy of the predicted 
small molecule hydration energies by the water models 
used in this study (in Tabic [TTJ) correlates well with their 
ability to reproduce experimental CHA. Namely, RMS 
errors against experiment in hydration free energy esti¬ 
mates for TIP5P, TIP3P, TIP4P-Ew and OPC are, re¬ 
spectively: 1.85, 1.10, 1.14, and 0.97 kcal/mol (based on 
20 neutral molecules, see Computational Methods, cho¬ 
sen [51] to span the entire hydration energy range of the 
set of 504 neutral molecules). The correlation suggests 
that the highly accurate experimental 77 * for the CHA- 
conjugate pairs of neutral small molecules can be used to 
calibrate and improve water models, especially their abil¬ 
ity to accurately describe microscopic electrostatic effects 
important for molecular hydration. The relative quantity 
77 * is expected to be insensitive to parameters describing 
the molecules in “CHA-conjugate” pair (such as the spe¬ 
cific parametrization of atomic partial charges) meaning 
that predicted 77 * mostly characterizes propensity of the 
water model to cause CHA. This propensity is strong, 
consistent with the large value of the charge-inversion 
symmetry breaking separation S, Fig. [2j comparable to 
the positive-negative charge separation in water models 
that approximate the experimental CHA reasonably well, 
Table [H] The CHA strength is also consistent with a rel¬ 
atively large value of the H 2 octupole moment of water 
molecule in liquid phase predicted by several independent 
quantum mechanical estimates [4611 5 21. 

As we have noted before, for each CHA-conjugate 
molecule pair one can extract AG sym directly from the 
experimental hydration free energies, which is essentially 
the average hydration energy of the pair. This quantity 
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can also be used to assess the quality of water models 
with respect to their ability to describe the symmetric 
response to solvated charge. For example, for the pairs 
in Table El OPC water model gives the best agreement 
A G sym cs -7.3 kcal/mol for each pairs, vs. the experi¬ 
mental value of ~ -8 kcal/mol; AG sym for other water 
models tested here are less accurate, Fig. 0] In the fu¬ 
ture, it would be interesting to test other water models, 
including polarizable ones, with respect to their ability to 
reproduce experimental 77* and A G sym . We expect that 
for the current polarizable models, for which the polariz¬ 
ability is added at the dipole level 0 , l53l - l55l ] , one would 
observe similar 77* to the corresponding non-polarizable 
“base” model to which the polarizability is added. This 
is because the dipole moment response is symmetric with 
respect to solute charge inversion. 



FIG. 6 : Charge distribution in OPC water model [SQL [5l(|. 
The distance between the oxygen center and the positive 
charges is noticeably shorter than the ~ 1 A assumed by 
common 3 point charge water models fid]. [42:]. 

Having quantified the intrinsic CHA effect of water by 
77* for a pair of neutral, cation-like and anion-like small 
molecules, we inquire if we can use the accurate exper¬ 
imental value of 77* to set a lower bound on the value 
of 77*“" for the (K + /F _ ) pair and thus reduce the dra¬ 
matic uncertainty associated with the original, ionic CHA 
shown in Fig. [T| To this end, we decompose the molec¬ 
ular 77* onto individual atomic contributions. To set the 
bound on 77* we single out the “CHA-dominant atom” 
for each molecule, Fig. [5J and use the fact that the CHA 
effect of these two atoms is larger than that of the com¬ 
bined effect from all the atoms in the molecular pair. 
This is because of partial cancellation of CHA for neu¬ 
tral molecules. Note that if the CHA-dominant atoms 
for these CHA-conjugate molecule pairs were exactly of 
the same size as those of the K + /F _ pair, then 77* from 
Table U| could be used directly to set the lower bound 
on the r]* lon . In reality, the radii of the CHA-dominant 
“cation” or “anion” atoms, Fig.[5j are different from that 
of K + or F _ , see SI. Nonetheless A A G scales with the 
ion size according to [2j which is consistent with experi¬ 
ment [23l425| | . We use this simple scaling relationship to 
relate atomic (CHA-dominant atoms of the molecules) 
and ionic 77*, see SI. To obtain more refined estimate of 
the bound, one needs to consider additional contributions 
specific to ionic AAG, specifically the non-polar interac¬ 
tions estimated in e.g. Ref. [56j and water-vapor inter¬ 
face potential. For the oppositely charged ions, both the 
liquid-vapor interface potential [30] and the ion cavity 


potential [57} contribute to AAG: the net contribution 
is | 2 e| <j>i n t, where <pi n t is the (positive) electrochemical 
potential corresponding to taking an ion from vacuum 
to bulk water. Note that the contribution is positive for 
the positive potential, and thus can only increase r]* lon 
relative to molecular rf. Experimental estimates of fint 
range from +25 mV [28} to +140 mV [29|. For the lower 
bound estimate of 77 *“", taking the lower bound on the 
experimental estimate of <f>i n t = 25 mV, results in a neg¬ 
ligible contribution to AAG. 

The end result is a lower bound of 77 * lOTl (K + /F _ ) set by 
the molecular rf: r]* lon > O. 5 O 77 *+ 0.14. Here, 0.14 is the 
non-polar contribution ton* 1 ™, which, according to a re¬ 
cent theoretical estimate[56} we use here, is not negligible 
for the K + /F _ ion pair, and originates mostly from the 
difference in ion-water dispersion interactions for these 
ions. Substituting 77 * = 0.48 from Table [H we arrive at: 

7 ?* ion ( K + / F -) > 0.38 (4) 

This lower bound estimate suggests that the “original” 
relative CHA effect for the K+/F - ion pair, Fig. [TJ is 
also likely to be strong. 

To conclude this section, we propose an experiment 
that can, in principle, estimate AAG and 77 * of K+/F - 
pair accurately, avoi ding the uncertainties discussed ear¬ 
lier, Fig. [T] The idea (58| is to measure hydration energies 
of two net neutral ion pairs, K + M+ and L + F _ , where L + 
and M+ are relatively large ions (compared to the K + , 
F _ ions) of similar size. For such large ions, |AG(L)| ~ 
|AG(M)| < |AG(K + )|,|AG(F“)|. Thus, in very dilute 
aqueous solutions of K + M - and L + F _ , AG(K + M~) — 
AG(L + F“) = AG(K + ) - AG(F") + (AG(M“) - 
AG(L+)) ~ AG(K+)-AG(F“) = AAG(K+/F“). Once 
AAG(K+/F“) is quantified in this manner, accurate 
rf = AAG(K+/F“)/(AG) can be obtained since the 
average (AG) for the ion pair can be quantified very ac¬ 
curately. 

Note that a variation on the above idea might be used 
to develop an alternative procedure for parameterization 
of ion force-fields for classical simulations f+l - ltiTT l . The 
first step would involve experimental determination of 
hydration energy of a salt nF - L n+ where L n+ is a large 
“composite” ion such that for its surface atoms the force- 
field parameters are well established and are not related 
to small ion parameters, e.g. Co(NH 3 )g 3 (cobalt hex- 
ammine) or Ni(NH 2 )g 2 ion. The next step would be to 
perform a series of explicit solvent simulations of nF _ 
L n+ neutral pair to fit F” parameters to match experi¬ 
mental hydration energy of 7 iF - L rl+ salt for each water 
model of interest. Once F _ parameters are known, one 
can perform a set of simulations to obtain cation (e.g. 
K + ) parameters by fitting against known ion pair hydra¬ 
tion energies. With the known cation parameters, the 
procedure can be continued to estimate the parameters 
for other anions. The oulined protocol would avoid a 
number of difficulties associated with both experimental 
and computational [38] estimates of hydration energies of 
charged species. 
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III. CONCLUSIONS 


In this article we have proposed an approach to use 
highly accurate experimental hydration energies of small 
neutral molecules to quantify the effect of charge hy¬ 
dration asymmetry (CHA) - the charge-asymmetric re¬ 
sponse of water to microscopic electric field. Character¬ 
istic dependence of hydration free energies on the sign 
of charged solutes of similar radii such as the K + /F“ 
pair is the best known manifestation of CHA. However, 
for ion pairs, the quantification of the hydration free en¬ 
ergy, and hence the associated CHA, is highly uncertain 
(~ 300% difference between four available comprehen¬ 
sive sets of experimental data) due to a variety of fun¬ 
damental, and technical difficulties [E], [2(1 Ho|, H&, 62j • 
Here we overcome these difficulties and accurately quan¬ 
tify CHA by proposing an approach that allows us to 
separate charge-asymmetric and charge-symmetric parts 
of the hydration free energies of neutral solutes. The 
crux of the approach is the identification of pairs of neu¬ 
tral solutes that show charge-asymmetric response of wa¬ 
ter similar to that that of the K + /F _ pair, based on 
a set of quantitative criteria we infer from the behav¬ 
ior of K + /F _ pair with respect to the charge-symmetry 
breaking perturbation. For this purpose we use a com¬ 
bination of “CHA-aware” implicit solvation model and 
free energy perturbation simulations in different explicit 
water models to search through a large, comprehensive 
set of small, neutral drug-like molecules. The search has 
yielded two “CHA-conjugate” pairs of neutral molecules 
that behave just like the K + /F^ pair with respect to 
asymmetric charge hydration. Unlike the correspond¬ 
ing quantity for K + /F _ pair, the measured relative hy¬ 
dration asymmetry of neutral solutes is very accurate. 
The difference between experimental hydration energies 
within each pair of these special neutral molecules, rela¬ 
tive to the pair’s average hydration energy, quantifies the 
intrinsic charge-asymmetric response of real water, free of 
extrinsic, complicating factors such as the energetic cost 
of charge species to cross the liquid/vapor boundary. 

Quantitatively, we find the asymmetry of the water re¬ 
sponse to hydrated microscopic charge (size about 1 A) 
to be close to one-half of the average charge hydration 
energy of the charge itself, which means that the charge- 
asymmetry of the intrinsic aqueous response to micro¬ 
scopic fields is strong. Given the paucity of available 
experimental characteristics of electrostatic properties of 
water molecule in liquid phase at microscopic level, this 
result is important. 

The availability of a novel, accurate reference value for 
the charge hydration asymmetry made possible through 
this work, should be of interest in its own right. For 
example, we have found that the ability of an explicit 
water model to predict the correct experimental relative 
CHA for just one pair of neutral solutes correlates well 
with the accuracy of the model in predicting absolute 
hydration energies of small neutral molecules covering a 
wide range of hydration energies. The observation sug¬ 


gests an immediate use of the proposed neutral molecule- 
based CHA reference for testing, and ultimately improv¬ 
ing, solvent models. Note that computational estimate 
of hydration energies of small neutral solutes is now in¬ 
expensive, straightforward and free from serious compli¬ 
cating issues associated with analogous computations for 
charged species. The inclusion of the proposed accurate 
CHA reference as an additional optimization target in 
the process of water model construction may eventually 
result in more accurate models of water. 


IV. COMPUTATIONAL METHODS 

A. The CHA-GB equation 

The implicit solvent hydration free energies were com¬ 
puted using the CHA aware generalized Born equa¬ 
tion [Til l 


A G=- 


1 

2 



QiQj 

rCHA—GB ’ 

Jij 


( 5 ) 


where ^ is the charge of atom i, f- 


CHA-GB 

ij 


RiR^e 


Ri is its effective Born radius, 


r'ij is the spatial separation between atom i and j, 
e = 80 is the dielectric constant of water. Here A, = 


Ri 1 + sgn 


E, 


Ri—Rs -\-pu 


introduces the 


CHA scaling of the effective Born radii, Ri analogous 
to the CHA scaling of the Born radii in [T] r is a posi¬ 
tive constant 0(1) that controls the effective range of the 
neighboring charges ( j ) affecting the CHA of atom (*). 
The effective Born radii Ri are computed using the nu¬ 
merical implementation of the “R6” formula [63j, over the 
dielectric surface (boundary). The later is obtained by 
increasing [jit the intrinsic atomic radii by R s = 0.52 A , 
see Ref [181] , and further rolling a probe of radius R w — R s , 
where R w = 1.4 A defines the size of water molecule. 


B. Selection of small molecules and identification 
CHA-conjugate pairs 

To identify CHA-conjugate molecule pairs from the 
original 504 molecule set j34| we first selected 250 most 
rigid molecules that undergo very little (<0.3 A) con¬ 
formational change, as seen from molecular dynamics 
trajectories [65j. We then picked 60 “polar dominant” 
molecules - those with the smallest nonpolar AG to po¬ 
lar AG ratio (<0.2), based on previous solvation energy 
estimates 0] in TIP3P water. We first performed single 
point AG estimates using a CHA-aware implicit solvation 
model; the polar part was modeled using the CHA-GB 
equation [5j and the nonpolar part modeled as 'jSASA , 








where SASA is the solvent accessible surface area com¬ 
puted numerically @ and 7 = 0.005 kcal/mol/A 1 2 3 4 5 6 7 8 . The 
topology and coordinates for the molecules were obtained 
from Ref. j34|. Intrinsic atomic radii set and r, [5j were 
optimized (using Nelder-Mead simplex algorithm H3) 
against the experimental solvation free energy for the 60 
polar molecules using different values of S in the range (0, 
1.0 A). For each value of <5 we performed 2000 separate 
optimizations using random initial guess; the parameter 
set corresponding to the smallest deviation from the ex¬ 
periments, was used for analysis. The resulting AG for 
different values of 5 were used to shortlist an initial se¬ 
lection of 11 molecule pairs (17 distinct molecules) that 
showed promising CHA-conjugate behavior with respect 
to the AAG(i5) gap. The initial set was then refined via 
careful explicit water free energy perturbation calcula¬ 
tions in TIP3P, TIP4P-Ew and TIP5P, resulting in two 
CHA-conjugate pairs, Table U that exhibit near perfect 
symmetric, monotonically increasing AAG(S) gap. 

C. Explicit solvation free energies 

Molecule topolo gy a nd coordinate files were prepared 
in an earlier work Tbit . using GAFF [68j] small molecule 
parameters assigned by Antechamber 14 |69j, and the 
partial charges were assigned using the Merck-Frosst im¬ 
plementation of AM1-BCC [70}. The hydration free en¬ 
ergy calculations in explicit water were performed in 
GROMACS 4.6.5 Q using standard free energy per¬ 
turbation (FEP) calculations 0] the coulomb and van 
der Waals coupling was reduced from 1 to 0 using 20 in¬ 
termediate A values. Molecules were solvated in triclinic 
box with at least 12 A from the solute to the nearest 
box edge. Real space electrostatic cutoff was 10 A, and 


long-range electrostatic interactions were calculated us¬ 
ing periodic boundary conditions via. the particle mesh 
Ewald (PME) summation [x2, 0| and all bonds were 
restrained using the LINCS algorithm. Production sim¬ 
ulations were 5 ns in length at each A value, and free 
energies and the associated uncertainties were computed 
using the Bennett acceptance ratio (BAR), namely the 
gbar feature in GROMACS 4.6.5. The FEP estimates 
were validated by ensuring convergence of forward (turn¬ 
ing the coupling on) and backward (turning the coupling 
off) computations for two molecules, and also by compar¬ 
ing with the FEP estimates for TIP3P in Ref. (34j. The 
average errors computed using gbar were roughly 0.05 
kcal/mol (see SI); the forward and backward computa¬ 
tions were comparable with a difference of 0.04 kcal/mol 
for the two molecules. The same topology and coordinate 
files were used along with identical simulation protocol 
to perform the FEP calculations in each of the four wa¬ 
ter models, TIP5P, TIP3P, TIP4P-Ew and OPC [Ho| for 
17 molecules (including the 11 shortlisted pairs, above). 
FEP calculations were also performed for 3 additional 
molecules, such that the set of these 20 molecules span 
the range of AG from -0.7 to -10.0 kcal/mol, which is 
close to the entire experimental range seen in the set of 
504 small neutral molecules. The estimated AG for the 
set of 20 molecules was used to compare the four wa¬ 
ter models against the experimental 77 *. Their estimated 
solvation free energies are provided in SI. 
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